home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / c / matrixc.zip / MATRIX.C < prev    next >
C/C++ Source or Header  |  1990-05-10  |  6KB  |  307 lines
  1. /* Matrix operation functions and typedefs */
  2. /* Written by Nigel Salt */
  3.  
  4. #include <stdio.h>
  5. #include <malloc.h>
  6.  
  7. typedef struct
  8. {
  9.   int rows;
  10.   int cols;
  11.   double *block;
  12. } matrix,*matrixptr;
  13.  
  14. /* Function prototypes */
  15. void mprint(matrixptr);
  16. void smmult(matrixptr,double);
  17. int madd(matrixptr m1,matrixptr m2,matrixptr dm);
  18. int mmult(matrixptr m1,matrixptr m2,matrixptr dm);
  19. int mcopy(matrixptr sm,matrixptr dm);
  20. int mtrans(matrixptr sm,matrixptr dm);
  21. double det(matrixptr m);
  22. int minv(matrixptr sm,matrixptr dm);
  23. int nsolve(int rows,double *data);
  24. int mid(matrixptr);
  25. void mzero(matrixptr);
  26.  
  27. /* function definitions */
  28. void mprint(m)
  29. matrixptr m;
  30. {
  31.   int i,j;
  32.   int cols;
  33.   cols=m->cols;
  34.   for (i=0;i<m->rows;i++)
  35.     {
  36.     printf("\n");
  37.     for (j=0;j<cols;j++)
  38.       {
  39.       printf("%8.2lf",*(m->block+i*cols+j));
  40.       }
  41.     }
  42. }
  43.  
  44. /* add two matrices giving a third matrix */
  45. int madd(m1,m2,dm)
  46. matrixptr m1,m2,dm;
  47. {
  48.   int i,j;
  49.   if (m1->rows!=m2->rows||m1->cols!=m2->cols)
  50.     {
  51.     fprintf(stderr,"\nmadd error - matrices different sizes");
  52.     return(1);
  53.     }
  54.   if (m1->rows!=dm->rows||m1->cols!=dm->cols)
  55.     {
  56.     fprintf(stderr,"\nmadd error - matrices different sizes");
  57.     return(1);
  58.     }
  59.   for (i=0;i<m1->rows;i++)
  60.     for (j=0;j<m1->cols;j++)
  61.       *(dm->block+i*m1->cols+j)=*(m1->block+i*m1->cols+j)+\
  62.       *(m2->block+i*m1->cols+j);
  63.   return(0);
  64. }
  65.  
  66. int mcopy(sm,dm)
  67. matrixptr sm,dm;
  68. {
  69.   int i,j;
  70.   if (dm->rows!=sm->rows||dm->cols!=sm->cols)
  71.     {
  72.     fprintf(stderr,"\nmcopy error - matrices different sizes");
  73.     return(1);
  74.     }
  75.   for (i=0;i<dm->rows;i++)
  76.     for (j=0;j<dm->cols;j++)
  77.       *(dm->block+i*dm->cols+j)=*(sm->block+i*dm->cols+j);
  78.   return(0);
  79. }
  80.  
  81. /* multiply matrix by scalar double */
  82. void smmult(m,s)
  83. matrixptr m;
  84. double s;
  85. {
  86.   int i,j;
  87.   for (i=0;i<m->rows;i++)
  88.     for (j=0;j<m->cols;j++)
  89.       *(m->block+i*m->cols+j)=*(m->block+i*m->cols+j)*s;
  90. }
  91.  
  92. int mmult(m1,m2,dm)
  93. matrixptr m1,m2,dm;
  94. {
  95.   int i,j,k;
  96.   double cellval;
  97.   if (m1->cols!=m2->rows)
  98.     {
  99.     fprintf(stderr,"\nmmult error - matrix 1 cols must = matrix 2 rows");
  100.     return(1);
  101.     }
  102.   if (m2->cols!=dm->cols)
  103.     {
  104.     fprintf(stderr,"\nmmult error - dest matrix cols must = matrix 2 cols");
  105.     return(1);
  106.     }
  107.   if (m1->rows!=dm->rows)
  108.     {
  109.     fprintf(stderr,"\nmmult error - dest matrix rows must = matrix 1 rows");
  110.     return(1);
  111.     }
  112.   for (i=0;i<m1->rows;i++)
  113.     for (j=0;j<m2->cols;j++)
  114.       {
  115.       cellval=0.0;
  116.       for (k=0;k<m1->cols;k++)
  117.         {
  118.         cellval+=*(m1->block+i*(m1->cols)+k) * (*(m2->block+k*(m2->cols)+j));
  119.         }
  120.       *(dm->block+i*dm->cols+j)=cellval;
  121.       }
  122.   return(0);
  123. }
  124.  
  125. int mtrans(sm,dm)
  126. matrixptr sm,dm;
  127. {
  128.   int i,j;
  129.   if (dm->rows!=sm->cols)
  130.     {
  131.     fprintf(stderr,"\nmtrans error - dest matrix rows must = source matrix cols");
  132.     return(1);
  133.     }
  134.   if (sm->rows!=dm->cols)
  135.     {
  136.     fprintf(stderr,"\nmtrans error - source matrix rows must = dest matrix cols");
  137.     return(1);
  138.     }
  139.   for (i=0;i<sm->rows;i++)
  140.     for (j=0;j<sm->cols;j++)
  141.       {
  142.       *(dm->block+j*dm->cols+i)=*(sm->block+i*sm->cols+j);
  143.       }
  144.   return(0);
  145. }
  146.  
  147. double det(m)
  148. matrixptr m;
  149. {
  150.   double p1,p2,p3,d;
  151.   int i,j,k;
  152.   if (m->cols!=m->rows)
  153.     {
  154.     fprintf(stderr,"\ndet error - matrix must be square");
  155.     return(0.0);
  156.     }
  157.   d=0;
  158.   for (i=0;i<m->cols;i++)
  159.     {
  160.     p1=p2=1.0;
  161.     p3=*(m->block+i);
  162.     k=i;
  163.     for (j=1;j<m->cols;j++)
  164.       {
  165.       k=(k+1)%m->cols;
  166.       p1*= *(m->block+j*m->cols+k);
  167.       p2*= *(m->block+(m->cols-j)*m->cols+k);
  168.       }
  169.     p3*=(p1-p2);
  170.     d+=p3;
  171.     }
  172.   return(d);
  173. }
  174.  
  175. int minv(sm,dm)
  176. matrixptr sm,dm;
  177. {
  178.   int i,j,k,l;
  179.     int nrow,ncol;
  180.   double *d;
  181.   if (sm->rows!=dm->rows||sm->cols!=dm->cols)
  182.     {
  183.     fprintf(stderr,"\nminv error - matrices must be same size");
  184.     return(1);
  185.     }
  186.   if (det(sm)==0.0)
  187.     {
  188.     fprintf(stderr,"\nminv error - matrix is singular");
  189.     return 1;
  190.     }
  191.   d=(double *)(malloc(sizeof(double)*sm->rows*(sm->cols+1)));
  192.   if (d==(double *)NULL)
  193.     {
  194.     fprintf(stderr,"\nminv error - insufficient memory");
  195.     return(1);
  196.     }
  197.   for (i=0;i<sm->rows;i++)
  198.     {
  199.         nrow=i-1;
  200.         ncol=i-1;
  201.     for (j=0;j<sm->rows;j++)
  202.       {
  203.       nrow=(nrow+1)%sm->rows;
  204.       if (j==0)
  205.         *(d+j*(sm->cols+1)+sm->cols)=1;
  206.       else
  207.                 *(d+j*(sm->cols+1)+sm->cols)=0;
  208.             for (k=0;k<sm->cols;k++)
  209.                 {
  210.                 ncol=(ncol+1)%sm->cols;
  211.                 *(d+j*(sm->cols+1)+k)=*(sm->block+nrow*sm->cols+ncol);
  212.                 }
  213.       }
  214.     if (nsolve(sm->rows,d))
  215.       {
  216.             fprintf(stderr,"\nminv error - cannot use nsolve on row %u",i);
  217.       free((char *)d);
  218.       return 1;
  219.       }
  220.     else
  221.             {
  222.             nrow=i-1;
  223.             for (j=0;j<sm->rows;j++)
  224.                 {
  225.                 nrow=(nrow+1)%sm->rows;
  226.                 *(dm->block+nrow*sm->cols+i)=*(d+j*(sm->cols+1)+sm->cols);
  227.                 }
  228.       }
  229.     }
  230.   free((char *)d);
  231.   return 0;
  232. }
  233.  
  234. int nsolve(rows,data)
  235. int rows;
  236. double *data;
  237. {
  238.   int i,j,k;
  239.   int cols;
  240.   double dtemp;
  241.   cols=rows+1;
  242.   for (i=0;i<rows;i++)
  243.     {
  244.     for (j=i;j<rows&&*(data+j*cols+j)==0.0;j++);
  245.     if (*(data+j*cols+j)==0.0)
  246.       {
  247.       fprintf(stderr,"\nnsolve error - singular matrix");
  248.       return 1;
  249.       }
  250.     if (j!=i)
  251.       {
  252.       for (k=0;k<cols;k++)
  253.         {
  254.         dtemp=*(data+i*cols+k);
  255.         *(data+i*cols+k)=*(data+j*cols+k);
  256.         *(data+j*cols+k)=dtemp;
  257.         }
  258.       }
  259.         for (j=cols-1;j>=0;j--)
  260.       {
  261.             *(data+i*cols+j) /= *(data+i*cols+i);
  262.       }
  263.     for (j=i+1;j<rows;j++)
  264.       {
  265.             for (k=cols-1;k>=i;k--)
  266.         *(data+j*cols+k)-=*(data+j*cols+i) * *(data+i*cols+k);
  267.       }
  268.         }
  269.     for (i=rows-2;i>=0;i--)
  270.         {
  271.         for (j=cols-2;j>i;j--)
  272.             {
  273.             *(data+i*cols+cols-1)-= \
  274.             *(data+i*cols+j) * *(data+j*cols+cols-1);
  275.             *(data+i*cols+j)=0;
  276.             }
  277.         }
  278.   return 0;
  279. }
  280.  
  281. int mid(m)
  282. matrixptr m;
  283. {
  284.   int i,j;
  285.   if (m->rows!=m->cols)
  286.     {
  287.     fprintf(stderr,"\nmid error - matrix must be square");
  288.     return 1;
  289.     }
  290.   for (i=0;i<m->rows;i++)
  291.     {
  292.     for (j=0;j<m->cols;j++)
  293.       *(m->block+i*m->cols+j)=0.0;
  294.     *(m->block+i*m->cols+i)=1.0;
  295.     }
  296.   return 0;
  297. }
  298.  
  299. void mzero(m)
  300. matrixptr m;
  301. {
  302.   int i,j;
  303.   for (i=0;i<m->rows;i++)
  304.     for (j=0;j<m->cols;j++)
  305.       *(m->block+i*m->cols+j)=0.0;
  306. }
  307.